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ABSTRACT 

Time delays between lensed multiple images have been known to provide an interesting probe of 
the Hubble constant, but such application is often limited by degeneracies with the shape of lens 
potentials. We propose a new statistical approach to examine the dependence of time delays on the 
complexity of lens potentials, such as higher-order perturbations, non-isothermality, and substruc- 
tures. Specifically, we introduce a reduced time delay of the dimensionless form, and explore its 
behavior analytically and numerically as a function of the image configuration that is characterized 
by the asymmetry and opening angle of the image pair. In particular we derive a realistic conditional 
probability distribution for a given image configuration from Monte-Carlo simulations. We find that 
the probability distribution is sensitive to the image configuration such that more symmetric and/or 
smaller opening angle image pairs are more easily affected by perturbations on the primary lens po- 
tential. On average time delays of double lenses are less scattered than those of quadruple lenses. 
Furthermore, the realistic conditional distribution allows a new statistical method to constrain the 
Hubble constant from observed time delays. Wc find that 16 published time delay quasars constrain 
the Hubble constant to be Hq = 70 ± 6kms~^Mpc~^, where the value and its error are estimated 
using jackknife resampling. Systematic errors coming from the heterogeneous nature of the quasar 
sample and the uncertainty of the input distribution of lens potentials can be larger than the sta- 
tistical error. After including rough estimates of the sizes of important systematic errors, we find 
Hq = 68 ± 6(stat.) ± 8(syst.) kms^^Mpc^^. The reasonable agreement of the value of the Hubble 
constant with other estimates indicates the usefulness of our new approach as a cosmological and 
astrophysical probe, particularly in the era of large-scale synoptic surveys. 

Subject headings: cosmology: theory — dark matter — distance scale — galaxies: elliptical and 
lenticular, cD — gravitational lensing 



1. INTRODUCTION 

It has been known that time delays between multiple 
images of strong gravitational lens systems offer an in- 
teresting method to measure the Hubble constant Hq^ 
the most fundamental cosmological parameter t hat gov- 
erns the length and time scale of our universe (IRefsdall 
119641) . A huge advantage of this method is that it 
does not rely on so-called distance ladder and can mea- 
sure the global Hubble constant independently of any 
local measurements. Motivated by this time delays 
have been meas ured in mor e than 10 lensed quasar sys- 
tems (see, e.g., iKochaneS [2006). The situation as it 
prese nts is, how e ver, so mewhat confusing and controver- 
sial. iKochaneS ()2002l . [2003 ) claimed from the analysis 
of several lens systems that the Hubble constant should 
be relatively low, Hq ^ 50 kms~^Mpc~^. The time 
delay of SDSS_J1650+4251 also prefers the low Hub- 
ble constant ([V uissoz et "aLll2007f) . On the other hand, 
iKoopmans et a l. (2003) performed systematic mass mod- 
eling of B1608-I-656 using all available data from radio 
to optical and found constrained the value of the Hubble 
constant to be Hq — 75lg kms^^Mpc^^. The analysis 
of the smallest separation lens B0218-I-357 yields Hq = 
78±3kms"iMpc^^ (.Wuck nitz et al.ll2004f ). By combin- 
ing tim e delays in 10 lensed quasar systems ISaha et al] 
(|2006D obtained Hq = 72j:^i kms-^Mpc-i. 

The large variation of derived values of the Hubble 
constant from time delays reflects the fact that time de- 
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lays are quite sensitive to mass distributions of lens ob- 
jects, which leads to strong degeneracies between lens 
mass profiles and the Hubble constant. The most fun- 
damental, mathematic ally rigorous dege neracy is the 
mass-sheet degeneracy (iFalco et al.lflQsHh : inserting an 
uniform sheet instead of decreasing the mass normal- 
ization of lenses changes estimated values of the Hub- 
ble constant while leaving unchanged the image observ- 
ables. This degeneracy implies that the derived values 
of the Hubble constant also degenerates with radial den- 
sity profiles of of lens galaxies (Rcfsdal & Surdei 199 j: 
Witt et al.1 1199 5: Keeton & Kochanck 1997; Witt et al. 
20001 ISahallFoOO,: .Tada fc Futamase ,2000: .Wucknitz 



2002t iTreu fc KoopmansI 120021: iKochanekl l2002l l2003l: 



Oguri & Kawano' 2003: Rusin et al." '2003"; 'Sch echten 
2005; Mortsell e t al. 2005; Kawano fc Oguri 2006|). In 
addition, it may d egenerate with the angular structure 
of le nses as well (|Zhao fc OinI l2003t ISaha fc Wilhamsl 
l2006l ). The sensitivity of time delays on mass profiles 
suggests that by assuming the Hubble constant, which 
can be determined by other methods (iFreedman et aP 
[2OOII: l^mTark et al.ll200fil : iSpergel et~l.ll2007D . we can 
put constraints on mass distribution s of lenses, m par- 
ticular radial density profiles fe.g.. lOguri et al.l 12004 
iKochanek et"ani2n06bl : Dobke fc Kingll2006D . 

One way to overcome the degeneracies is to adopt lens 
systems whose lens potentials can well be constrained by 
many observational constraints. An example of such con- 
straints comes from lensed images of quasar host galaxies. 
The lensed host galaxies often form Einstein rings, which 



2 



Oguri 



accurately and independently determines the structure of 
the lens pote ntial as well as the shape of the lensed host 
galaxy (e.g.. iKeeton et al.l l200d: iKochanek eFall 120011 : 
iKoopmans et al.1 l2003f) . ^The small-scale structure of 
lensed quasars, such as radio jets and sub-components, 
determine lens potentials in a similar way as E instein 
rings (e.g.. lCohn et al.ll200lHWucknitz et al.ll2004f ). Fur- 
thermore, the measurement of velocity dispersions of 
lens galaxies serves as useful constraints on mass dis- 
tributions, hel ping to b reak degeneracies between mass 
models fe.g.. Treu fc Ko opmans 2002). Lens systems for 
which these strong additional constraints are available, 
sometime referred as "golden lenses" , have been thought 
to be an effective probe of the Hubble constant. 

Another potentially powerful, but less studied, method 
to measure the Hubble constant is a statistical approach. 
Even if each individual lens lacks strong constraints that 
allow detailed investigations of the mass distribution, by 
combining many lens systems one can put tight con- 
straints on the Hubble constant. As mentioned above, 
thi s approach was in some sense demonstrated recently 
by ISaha et all (|2006[ ) who combined 10 lensed quasar 
systems to constrain the Hubble constant. A caveat is 
that lensed quasars sometimes suffer from selection ef- 
fects. For instance, brighter lensed quasars with larger 
image separations are more likely to lie i n dense en- 
vironments such as groups and clusters (jOguri et al.l 
I2005t IOgurill2006[ ). thus the Hubble constant from those 
lensed quasars may be systematically higher than the 
true value without any correction of the effect of dark 
matter along the line of sight. This indicates the impor- 
tance of well-defined statistical sample of lensed quasars 
for which we can quantitatively estimate and correct 
the selection effect. While the statistical sample has 
been very limited so far, containing ^ 20 lenses e ven for 
the la rgest lens sample (Mvcrs ct al. 2003; Brown e" et al.l 
I2003D . larger lens surveys will be soon available by ongo- 
ing/future le ns surveys such as done by the Sloan Digital 
Sky Survey (jOguri et al.ll2006l ). the Large Synoptic Sur- 
vey Telescope (LSST)^, and the Supernova/ Acceleration 
Probe (SNAP) . Future lens surveys will als o find strong 
lensing of distant supernovae (e.g., iQguri et all 120031) 
for which time delays can easily be measured (but see 
iDobler fc Keetonl l2006l ). Therefore the statistical ap- 
proach is growing its importance. 

In this paper, we study how time delays depend on 
various properties of the lens potential and image config- 
urations, which is essential for the determination of the 
Hubble constant from time delays. Using both analytic 
and numerical methods, we show how time delays are af- 
fected by several complexity of the lens potential, such as 
radial mass profiles, external perturbations, higher order 
multipoles, and substructures. We then derive the ex- 
pected distributions of time delays by adopting realistic 
lens potentials. The distribution, in turn, can be used to 
derive statistically the value of the Hubble constant from 
observed time delays . This appr oach differs from the sta- 
tistical argument bv lSaha et al.l (2006) in the sense that 
they fit image positions of individual lens systems to con- 
strain the Hubble constant and then combined results of 
all the lens systems: Our new approach does not even 

^ http://www. lsst.org/ 
^ http://snap.lbl.gov/ 



require modeling of each lens system. This has an advan- 
tage that we can include lens systems that have too few 
constraints to determine the lens pot ential. In thi s sense 
the approach is extension of study bv lOguri et al.l ()2002f ) 
in which only spherical halos are considered to compute 
time delay distributions. We note that the metho dology 
is similar to that adopted bv lKeeton eFall (|2003l l2005h 
who derived distributions of flux ratios of image pairs to 
identify small-scale structure in lens galaxies. 

In addition to the measurement of the Hubble con- 
stant, the sensitivity of time delays on mass models, 
which we explore in this paper, offers guidance on the 
usefulness of each time delay measurement as a cosmo- 
logical or astrophysical probe: If time delays at some 
image configuration is quite sensitive to detailed struc- 
ture of the lens such as higher-order multipole terms with 
small amplitudes or substructures, which are difficult to 
be constrained even for best-studied lens systems, it is al- 
most hopeless to use these time delays to extract either 
radial mass profiles or the Hubble constant. Our result 
can be used to assess quantitatively which lens systems 
are more suitable for detailed studies, i.e., less sensitive 
to the complexity of the lens potential . There has been 
several insightfu l analytic work fe.g.. I Witt et al. I [20001 : 
iKochaneS |2002| ) , but here we perform more systematic 
and comprehensive survey of model dependences of time 
delays by parameterizing image configurations of lensed 
quasar systems using dimensionless quantities. 

This paper is organized as follows. In S|2] we intro- 
duce several dimensionless quantities that are used to 
explore the model dependence of time delays. We study 
the sensitivity of time delays on the various lens poten- 
tials analytically in ^J3l Section [4] is devoted to construct 
the conditional distribution of time delays from realistic 
Monte-Carlo simulations. We compare it with observed 
time delays in f|5l and constrain the Hubble constant 
in fjni Finally discussion of our results and conclusion 
are given in SJT] Throughout the paper we adopt a flat 
universe with the matter density = 0.24 and th e 
cosmological constant — 0.76 (jTegmark et al.|[2006[ ). 
although our results are only weakly dependent on the 
specific choice of the cosmological parameters. The Hub- 
ble constant is sometimes described by the dimensionless 
form h = iJo/(100kms-iMpc^i). 

2. CHARACTERIZING TIME DELAY QUASARS 

Let us consider a lens system in which a source at u = 
(u, v) is multiply imaged at the image positions = 
{xi, Hi). We also use the polar coordinates for the image 
positions, = {xi,yi) = (r^ cos 6*^, sin 6*^). We always 
choose the center of the lens object as the origin of the 
coordinates. Time delays between these multip ly images 
are given by fe.g. JBlandford fc NaravMlll986f ) 

. _ 1 + DqiDqs 
"'^ 2c As 
X [(x, - u)2 - (xj - u)2 - 20(xO + 20(x,)] , (1) 

where z; is the redshift of the lens, c is the speed of light, 
and Do\, D^s, and Dig are angular diameter distances 
from the observer to the lens, from the observer to the 
source, and from the lens to the source, respectively. The 
lens potential 0(x) is related to the surface mass density 
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of the lens S(x) by the Poisson equation: 



(2) 



with Scrit = c^Dos/{4:TrGDo\D\s) being the critical sur- 
face mass density (G is the gravitational constant). Note 
that image positions and source positions are related by 
the lens equation 



u 



(3) 



Equation ([T]) involves unobservable quantities such as 
u and (/'(x), indicating that time delay s in general d epend 
on details of mass models. However, IWitt et all ()2000D 
has shown that for generalized isothermal potential 



(t>{x.) = rF{0), 



(4) 



where F(9) is an arbitrary function of 9, time delays can 
be expressed in a very simple form involving only the 
observed image positions: 



1 + zi DoiDo 



2c 



Du 



(5) 



where denote the distance of image i from the center of 
the lens galaxy. Motivated by this analytic result, in this 
paper we consider the reduced time delay that is defined 
by: 



2c 



1 + zi DoiDo 



(X, - U)2 - (Xj - U)2 - 2(l>{Xi) + 2(l>{Xj) 



(6) 



In this definition, the reduced time delay S is always 
unity if the lens potential can be expressed by equa- 
tion ([?]), but can deviate from unity if the lens potential 
has more c omplicated struct ures. In particular, from the 
analysis of lKochanetJ ()2002f ) we obtain S — 2(1 — (k)) at 
the lowest order of expansion, where (k) is the average 
surface density in the annulus bounded by the images. 
This indicates that the deviation from the isothermal 
mass distribution has a direct impact on the reduced time 
delay. But as we will see later it is affected to some extent 
by other factors such as external perturbations or small- 
scale structures as well. Thus we can regard the reduced 
time delay as a measure of the complexity of the lens. 
In addition, equation ^ indicates that we can compute 
reduced time delays for observed lensed quasar systems 
with measured time delays by assuming the value of the 
Hubble constant. In this sense, the reduced time delay S 
is a key quantity that represents a link between measured 
time delays of lens systems and theoretical lens models. 
We point out that 5 is dimensionless because time de- 
lays are proportional to the square of the size of a lens 
(the Einstein radius). This allows us to directly compare 
values of reduced time delays for different lens systems 
that have differ ent sizes. 

We note that ISahal (|2004f ) adopted similar but differ- 
ent dimensionless quantity to explore the dependence of 
time delays on mass models. In the paper, the parame- 
ter essentially same as equation ^ was also considered, 
but it was discarded because of the correlation with time 
delays. However, the correlation just reflects the effect 



of surrounding dark matter that is larger for wider sep- 
aration lenses (jOguri et al.l 120051 : IOgurill2006f ). Put an- 
other way, such correlation is naturally expected from 
very different mass distributions and environments be- 
tween small (~ 1") and large (> 3") separation lenses. 
Indeed the apparent lack of the corre lation betwe en time 
delays and the scaled time delays in lSahal (|2004f) comes 
mostly from the large scatter among different image con- 
figurations: The effect of external mass is hindered by 
the large scatter his parametrization involves. Therefore 
in this paper we propose equation ([6]) as useful quantity 
to extract the mass dependences on time delays. 

Previous analytic calculations of time delays sug- 
gest that the model dependence of time delays is a 
strong function of image configurations (jWitt et al.ll2000l : 
Koclianck 2002). In this paper we characterize image 
configurations by the following two parameters. One is 
the asymmetry of the images define by 



R,j = 



(7) 



Again, Rij is dimensionless and does not depend on the 
size of the lens: Rij ~ means the images are roughly 
at the same distance from the lens galaxy, while Rij 1 
indicates very asymmetric configurations that one image 
lies very close to the lens center and the other image is 
far apart from the lens. The other parameter we use is 
the opening angle of images: 



(8) 



In this definition, if the images are directly opposite each 
other, 9ij ^ 180°. On the other hand, close image pairs 
such as merging images near cusp and fold catastrophe 
have 9ij ~ 0° . Note that both Rij and 9ij are observables 
in the sense that they can be derived without ambiguity 
for each observed lens system as long as the lens galaxy is 
identified: We do not have to perform mass modeling to 
compute these quantity from observations. In summary, 
our task of this paper is to explore model dependences 
of reduced time delay S as a function of image configu- 
rations parameterized by Rij and 9ij . 

3. ANALYTIC EXAMINATION 

In this section, we analytically examine the behavior 
of the reduced time delay S (eq. [6]) for various lens po- 
tentials, before showing the distribution of S for realistic 
complicated mass distributions. For this purpose, it is 
convenient to study in terms of multipole expansion: We 
consider the lens potential of the following from 



(3 



0, 



(9) 



where c„ is the dimensionless amplitude and 0„ is the po- 
sition angle. The coefficients are chosen such that i?Ein 
becomes the Einstein radius of the system if the am- 
plitude of the monopole term is cp = 1. Note that an 
external shear (e.g., Kochanek. 19911: iKeeton et al.lll997b 
can be described by this expression as /? = n = 2 and 
c„ = —7. For this potent ial the reduced time delay is 
given by (|Witt et al.ll2000l ): 



1 + ^2(1-/3)- 



(Xj) - 0(Xj) 



(10) 
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This can be rewritten as 



2-REii 



2-/3 



where S and X{Rij, 9ij) are defined by 

r: 



X{Rij , 0ij) cos{n0n - S) , 
(11) 



tan (5 i 



sin n9j 



sin n^,- 



cos n^j- 



rf cos nOi 



(12) 



-2(l-i?^^-)'^cosn6'y] 



1/2 



(13) 

Here we assumed Tj > without loss of generafity. Note 
that Rij and 9ij were defined in equations ([7]) and ([5]), 
respectively. 

The above expression of S (eq. [TT]) has several im- 
portant implications. First, On comes only in the last 
cosine term, thus assuming n 7^ and c„ is small S > 1 
and S < 1 occurs equally if 9n is uncorrelated with the 
image configuration (as we will see later, however, this 
is not necessarily the case). This indicates that we can 
reduce the effect of multipole terms by averaging many 
lens systems, suggesting the usefulness of our statistical 
approach. Second, since 2R^/{rj -\-ri) ^ 1 in most cases, 
the dependence of the potential {(3) and image configu- 
rations {Rij and 9ij) on the deviation from S = 1 is 
encapsulated in X{Rij,9ij), aside from the overall am- 
plitude c„. In what follows, we use X{Rij,6ij) to study 
the behavior of S. 

First, we consider quite symmetric cases {Rij — > 0) 
for which images are nearly the same distance from the 
lens center. Depending on the opening angles, limiting 
behaviors are given by 



X{Rij 0, 6ij) 



1-/3 
(3 Rij 



(cos nOij 
(cos nOii 



1) 
0) 



(14) 



{cos n6ij = — 1) 



Therefore, X{Rij, Oij) diverges unless cos nOij = 1. More 



rigorously, we should take the limit of i?^ 
0: 



and 9i. 



X{R, 



0) « (1 - /3) 



1 + 



1/2 



indicating that the divergence at the symmetric limit 
can be avoided if Oij ^ Rij, or more appropriately 
1 — cos nOij <C R^j. Since close image pairs are always 
near the circumference of the critical curve that has a 
nearly circular shape centered on the lens galaxy, such 
pairs in general have Rij <^ 9ij, indicating the diver- 
gence cannot be avoided for two close images. On the 
other hand, the divergence may not occur for opposite 
images 



180°), but only if n is even number. 
Inversely, if images are very asymmetric {Rij 
X {Rij, Oij) reduces to the following simple form 



X{Rij 1, Oij) 



2'3-i(l-/3) 



1), 



(16) 
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Fig. 1. — The behaviors of X{Rij ,9ij) (oq. [13]) as a function of 
the asymmetry parameter Rij. From thick to thin lines, the radial 
slope /3 is changed from 0.5 to 4.0. From top to bottom panels, the 
opening angle 8ij is fixed to cosnOij = —1, 0, and 1, respectively. 

This docs not depends on Oij, thus the opening angle is 
no longer important in this situation. However it shows 
stronger dependence on the radial slope (3. 

Finally we plot X {Rij, Oij) for various parameter val- 
ues in Figures [1] and [21 In Figure [TJ X {Rij, Oij) is plot- 
ted as a function of the asymmetry Rij . It grows quite 
rapidly as Rij approaches to zero for cos nOij ^ 1 . It 
also shows stronger sensitivity on the slope (3 at larger 
Rij . As is clear in Figure [51 the opening angle Oij be- 
comes more important for more symmetric lenses. All 
these behaviors are consistent with analytical arguments 
presented above. 

4. TIME DELAYS IN REALISTIC LENS MODELS 

In this section, we consider more realistic situations to 
in order to study expected spread of the reduced time de- 
lay S as a function of image configurations. Specifically 
we adopt theoretically and observationally determined 
distributions of lens potentials such as ellipticities, ex- 
ternal shear, substructures, and multipole components 
to make predictions on realistic probability distributions 
of S. The m ethodology is similar to that in lKeeton et al.l 
((200ai2005h who studied cusp and fold relations to iden- 
tify lenses with small-scale structure. 

4.1. Input Models 
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Fig. 2. — Same as Figure [T] but X(Rij,9ij) is plotted as a 
function of the opening angle 9ij while fixing Rij = 0.1 {top panel), 
0.5 {middle panel), and 0.9 {bottom panel). 

As primary lens galaxies we only consider elliptical 
galaxies because most (> 80%) of quas ar lenses are 
caused by massive elliptical galaxies (e. g., iTurner et al.l 
[T98llKoclia^[200llMoller et al.l[200l . We model the 
lens galaxy as a power-law elliptical mass distribution. 
The surface mass distribution is given by 



Vl -e cos 2(61 -6ie) 



2~a 



(17) 



where i?Ein is the Einstein ring radius and 9e is the po- 
sition angle of ellipse. The case a = 1 corresponds to 
the standard singular isothermal mass distribution. The 
ellipticity e, which is defined by e = 1 — g with q being 
the axis ratio of the ellipse, is related to e by 



l-(l-e)2 



l + (l-e)2 

The corresponding lens potential can be described by 

1 



(18) 



(19) 



with G{9) being the complex function of 9, but notice 
that G{e) = 1 if e = 0. 

Many previous work has shown that lens galax- 
ies indeed have nearly isothermal mass distribution. 



iRusin fc Kochanekl (|2005D obtained a = 0.94 ± 0.17 by 
combining the Einstein radii of many lensed quasars with 
the fundamental plane relation of elliptical galaxies. The 
slope of the galaxy density profile a — 1.09 ± 0.01 con- 
strained from the faint third image of PMN J1632— 0033 
is co nsistent with the nearly isothermal density pro- 
file dw mn et al.l I2004D . Detailed mass modeling of 
B1933-I-503 also indicate s nearly isothermal profile of the 
lens galaxy (jCohn et al.| l2001). By using measured veloc- 
ity dispersions of le ns galaxies of s everal lensed quasars , 
iTreu k, KoopmansI (|2004.. see also iHamana et a"l]|2007D 
derived a = 1.25 ± 0.2. iKoopmans et all (|2006f ) put 
tighter constrains a = 0.99^g'Q2, but the results are 
derived from a sample of much lower redshift lens sys- 
tems than typical time delay quasars. From these re- 
sults, in this paper we adopt the Gaussian distribution 
a = 1 ± 0.15 as a conservative input distribution of the 
slope. For the ellipticity, we use the Gaussian distri- 
bution with median e — 0.3 and dispersion 0.16 that 



contour shaoes of elliotical ealaxies (iBender et al.l 


1989 


Saglia et al. 


1993 


JJorgensen et al.lll995l;lRest et al.l 


2001 


Sheth et all 


2003 


). 



To allow more complex mass distribution of the galaxy, 
we add higher order multipolc terms to the potential: 



6m (x) = ii?|7 



r^'^il - rri^) Am cos m{e -9 m)- (20) 



The factor 1 — m? is inserted such that Am denotes 
the standard parametrization for the deviation of the 
mass density from an ellipsoid. We include only m = 3 
and 4 terms, because m > 5 perturbations have gen- 
erally not been reported. For both A^ and ^4, the 
amplitudes are distributed by the Gaussian with mean 
zero and dispersion 0.0 1 that is roughly co r isistent with 
reported distributions (jBender et al.1 119891 : ISaglia et al.l 
119931 iRest et al.ll2001[ ). It is also compatible with the 
level of the deviation i nferred from in dividual modeling of 
lensed quasar systems (| T rotter et aL|[20 00: Kawano et al 



2004; Kochanek & Dalar'200i ICongdon fc Keetonll200 
Yoo et al.. 2006) . We assume that position angles of these 



multipole perturbations and that of an ellipsoid are un- 
correlated with each ot her. In addition, we neglect the 
correlation of e and A^ (iKeeton et al.l 120031 ) for simplic- 
ity. Note that these multipole perturbations do not affect 
S if the lens galaxy is isothermal, a = 1, but can change 
quantitative results when combined with non-isothermal 
profiles and/or other small-scale structures. Since the 
effect of Ai itself is small, we expect the effect of the 
correlation between e and ^4, which we have neglected 
here, is also small. 

We also include external perturbations that are known 
to be important for individual mass modeling (e.g., 
iKeeton et al.lll997[ ). The potential of lowest order ex- 

^ The adopted d i stribu tion does not include the mean value of 
ITreu &: KoopmansI I I2004I) within its l-cr uncerta inty. However, we 
note that the result of lTreu &: KoopmansI 1 )20041) were drawn from 
only five lensed quasar systems and therefore scatter should be 
large. Moreover, many of lens systems used by ITreu &: KoopmansI 



II2004I ) appear to reside in dense environments, which may explain 
their large value of a. Indeed, Rusin & Kochanek (2005) used a 
larger sample of 22 lenses to derive the slope that is more consistent 
with our input distribution. 
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Fig. 3. — Dependences of the reduced time delay S on several lens potentials as a function of asymmetry Rij {left panels) or opening angle 
6ij (right panels). In each panel, reduced time delays of 500 double lenses (filled squares) and 500 quadruple lenses (open circles) obtained 
by the Monte-Carlo simulat ions are plotted. Magnification bias is not consid ere d at this stage. Fro m t op to bottom panels, we consider 
the external shear (/)e2 (eq. |21j'). th ird order external perturbation <f>^3 (eq. \22\). subhalos <f>g (eq. |24j ) . and non-isothermality a ^ 1 in 
the primary lens potential <f>Q (eq. 1191 ). Besides the bottom panels, we adopt an isothermal elliptical lens as primary lens galaxies (eq. 
|19) . a = 1) and ignore the other perturbations (e.g., in the second row we only consider 0e3 and ignore the external shear and subhalos), 
thus the effect of each potential can be measured by the deviation from H = 1. In all simulations higher order multipoles (eq. 1201 ) are 
included. 



ternal shear is given by 

0E2(x) = - 



cos 2(61 



(21) 



We adopt a log-normal distribution with median 7 = 
0.05 and dispersion 0.2 dex for the distribution of shear 
amplitude. It is consistent with expected shear distribu- 
tion from A^-body sim ulations ( Holder fc Schec htcr 20l)3'; 
iDalal &: WatsonI 12004) . In addition to exte rnal shear, 
we consider third order external perturbation ijKochanekl 
119911: iBernstein fc Fischal 1 1999( 1 



'/>E3(x) 



:Lr^ \cos( 



9^) - cos3{9 - 6,)] . (22) 



Here we assumed the external pcrturber is a singular 
isothermal object to relate the amplitudes and position 
angles of n = 1 and n = 3 terms. In general we have 
cr « 7^ ([Bernstein fc Fischeij Il999f ). thus for the am- 
plitude we adopt a log-normal distribution with median 



(7 = 7 and dispersion 0.2 dex. The adopted amplitude 
is roughly consistent with values obtained from mass 
mode ling of individual lens systems (e.g., Kawano et al.l 
120041 ). We allow small mis-alignment of the position an- 
gle 9a-, which could happen when the external perturber 
is not spherical, by adopting the Gaussian distribution 
around 9-y with dispersion 10°. 

Finally we consider substructures in the lens galaxy. 
A significant fraction of substructures (subhalos) in 
galaxies have been predicted as a na tural conse- 
quen c e of cold dark m atter cosmology (jMoore et al.l 
Il999t iKlvpin et al.l Il999f ). Anomalous flux ratios ob- 
served in many gravitational lens s ystems indicate 



that substrtictures are indeed present (Mctcalf fc Madau 



2001 



IChiball2002l: IDalal fc Kochan ek''2002 : Bradac et al 



2002i iKeeton et al l 12005 



.Ko chanck fc Dalall l20'ol 

Metcalf et al.ll2004l : IChiba et al.ll2005i) . Note that small 
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Fig. 4. — Contours of constant conditional probability p(S|/?ij) (left panels) and p(S|6ij) {right panels) computed from Monte-Carlo 
simulations of realistic lens potentials. The contours are drawn at the 68% (solid lines), 95% (dashed lines), and 99.7% (dotted lines) 
confidence levels. Upper panels show the conditional probability for double lenses, whereas lower panels are for quadruple lenses. For 
double lenses we show p(S|Sij) only at 8ij > 110° because we have too few image pairs with 6ij < 110° to construct the probability 
distribution. 

perturbations may also come from small halos along the to bk through bk/Rs oc (v/V)'^. We distribute subhalos 

randomly with an uniform spatial density in the pro- 
jected two-dimensional plane in order to take account of 
the s uggested anti-bias of the subhalo spatial distribu- 
tion (iDe Lucia et al.ll200l iMao et al.|[200S iQguri fc Lei 



line-of-sight i 


Keetonll2003l:lChen et al.ll2003HOEurill2005l: 


iMetcalfl 12005 


). Although time delays are thought to be 



less sensitive to small-scale structure than flux ratios that 
are determined by the second derivative of the time delay 
surface, they might be affected to some extent particu- 
larly when two images are close to each other. We model 
each subhalo by pseudo-Jaffe (truncated singular isother- 
mal) profile. The lens potential of this profile is 



</'PJ,fc(x) =&/c 



+ ak In r 



120041) . The resulting mass fraction of subhalos at around 
the Einstein radius is ~ 0.5%, being consistent with the 
expectation fr om iV-body simulations and analytic cal- 
culation (e.g.. IMao et"alll2004l : IOgurill2005f) . The effect 
of subhalos is described by the sum of the lens potential 
of each subhalo: 

0s(x) = <?!)pj,fc(x - Xsub,fe) - l:r^Ksuh- (24) 



,(23) 



fc 



where Uk is a truncation radius and bk is a mass nor- 
malization that coincides the Einstein radius for suf- 
ficiently large Uk- We adopt Uk = VbkRE assuming 
the truncation radius of the tidal ra dius of the sub- 
halo (see, e.g.. iMetcalf fc Madaull200lD . For the veloc- 
ity distribution we assume N{> v) = {lOv /V)~'^''^ in- 
side three times the Einstein radius of the lens galaxy, 
where v and V are velocity dispersions of the subhalo 
and halo, respectively. The velocities can be converted 



Here we subtracted the convergence averaged over all 
subhalos, Ksub, to conserve the total mass and radial pro- 
file of the lens galaxy. We include only 100 most massive 
subhalos mainly for computational reason, but small sub- 
halos are expected to have only very small effect on time 
delays. 

4.2. Simulation Method 

Using the model described above, we perform a large 
Monte-Carlo simulation containing > 10^ lensed image 
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Fig. 5. — Contour plot of median (left panels) and scatter {right panels) of the conditional probability p{S[Rij, 6ij) in the Rij-dij plane. 
Here the scatter is defined by the 68% confidence interval width in units of logS. The probability distribution for double {upper panels) 
and quadruple {lower panels) lenses are shown separately. Thick solid lines indicate the limit beyond which the number of image pairs in 
Monte-Carlo simulations is too small to construct the conditional probability. 



pairs. The simulation is done by first generating a lens 
potential according to the distributions summarized in 
t j4.1l All lengths are scaled by the Einstein ring radius 
i?E: Since our results do not depend on adopted length 
scales, we fix i?E to unity in the simulations. We gener- 
ate 10000 different lens potentials in total. For each lens 
potential, we place random sources with an uniform den- 
sity of ~ 100i?£^ i n the source p lane. We use a public 
software lensmodel (|Keetonll2001| ) to solve the lens equa- 
tion and compute time delays between multiple images. 
The uniform sampling in the source plane indicates that 
each lens potential is automatically weighted by the lens- 
ing cross section (see iKeeton et al.. 20031. To account for 
magnification bias as well, for each source we compute 
the total magnification factor /xtot, and when construct- 
ing probability distributions of reduced time delays (see 
below) we include a weight of /i^^^ , where g is a power- 
law slope of the luminosity function of source quasars. 
We adopt q = 2.1 that is releva nt for lenses identified by 
the CLASS (|Mvers et al.l 120031) . 



From the ensemble of image pairs, we compute condi- 
tional probability distribution functions of the reduced 
time delay S for given values of the asymmetry Rij 
and/or opening angle 9ij, i.e., p(S|i?jj), p{B\6ij), and 
p(S I Rij , 9ij ) . The probability distributions are computed 
separately for double and quadruple lenses to see how 
different distributions they exhibit. In what follow we 
ignore central faint images that are unobserved in most 
cases. 

4.3. Contribution of Each Potential 

Before presenting results that include all potentials, it 
is useful to see how each potential affects the distribu- 
tion of S. To see this we adopt isothermal elliptical lens 
potential plus multipole terms (0g + (t>M, « = 1), and 
add each lens potential: Since the isothermal galaxy al- 
ways yields S = 1 (see Sj2|) the effect of each potential 
term can be estimated by the deviation from S = 1. We 
also study the effect of non-isothermality {a ^ 1). The 
results are summarized in Figure [3] First, the effect of 
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Fig. 6. — Same as Figure |4] but observed values of S with error-bars (see Table [TJ are plotted on the contours by filled triangles with 
error-bars. The Hubble constant of h = 0.73 is assumed when computing H from observed time delays. 



external shear appears as a scatter around ^ = 1, which 
is consistent with discussion above O as well as that of 
iWitt et al.l (|2000D . Both the diverging scatter at R^j 
and the modest increase of the scatter at smaller 9ij are 
expected from our analytic examination. The scatter in- 
troduced by the third order external perturbation and 
subhalos is smaller than that of external shear, but it is 
still noticeable particularly for symmetric configurations. 
This means that even if we can estimate external shear 
and its orientation accurately for individual lens system 
by detailed mass modeling, these third order perturba- 
tions or subhalos can shift time delays, although these 
small perturbations mainly affect image pairs with small 
opening angles and only increase the scatter rather than 
shift the mean. We comment that the use of image pairs 
with small opening angles is also limited by the larger 
uncertainties of their time delay measurements. As ex- 
pected, non-isothermality has a large impact on S, but 
the size of its scatter is less dependent on image config- 
urations compared with other potentials. 

It is worth noting that there are systematic deviations 
from S = 1 for a few specific image configurations. For 
instance, at small 9ij external shear preferentially pro- 
duces image pairs with time delays smaller than the 
isothermal case. In addition, very asymmetric image 
pairs (Rij ^ 1) tend to have larger time delays when we 



consider non-isothermal lens galaxies. Such systematic 
shifts were not predicted in our analytic arguments in Sj3l 
thus invite careful consideration. We find that the selec- 
tion effect can explain these systematic shifts. For very 
asymmetric pairs, such configurations are possible only 
when the lens galaxy has a profile steeper than isother- 
mal for which inner critical curve is degenerated at the 
center. If the density profile is shallower than isothermal, 
any images (expect central faint images that we ignore 
in this paper) are outside the inner critical curve. This 
means that very asymmetric configuration, in which one 
image lies very close to the center of the lens potential, 
never occurs for the profile shallower than isothermal. 
Since shallower profiles correspond to smaller S, this in- 
creases the mean S at i?^ ~ 1. The effect of external 
shear is more complicated, but can be explained as fol- 
lows. Consider a close pair of two images. If we add 
external shear whose position angle direction is perpen- 
dicular to the segment that connects the two images, it 
separates two images. On the other hand, external shear 
that is parallel to the segment makes the two images 
closer. Therefore, combined with the steep dependence 
of the frequency of close image pairs on the opening an- 
gle, this results in the situation that at fixed small open- 
ing angle the direction of external shear is more likely 
to be parallel to the segment than perpendicular. The 
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discussion in S|3] indicates that the paraUel shear yields 
smaller time delays S < 1 for fixed image positions, thus 
we can conclude that close image pairs have statistically 
smaller time delays than those expected without external 
shear. 

4.4. Conditional Probability Distributions 

We are now ready to compute conditional probabil- 
ity distributions of the reduced time delay S with all 
lens potentials presented in ^ 34. II Figure [4] shows con- 
tours of constant conditional probability projected to 
one-dimensional surface, p{E\Rij) and p{E\9ij). The be- 
haviors of contours can well be understood from discus- 
sions in S|3] and ^ ^4.31 In particular, systematically small 
S at small opening angles and large S for very asymmet- 
ric image pairs, which we discussed in M.S) are clearly 
seen. It appears quadruple lenses have larger scatter 
on average: This is because strong perturbations (exter- 
nal shear, third order perturbation and subhalos), which 
cause the scatter around S = 1, also enhance the lensing 
cross section for quadruple images. 

To study the conditional probability distribution as a 
function of both Rij and 0^ , p{E\Rij ,9ij), in Figure [5] 
wc plot contours of median and 68% confidence inter- 
val of the probability distribution in the Rij-dij plane. 
The features we have discussed, such as the divergence 
of scatter at Rij — + 0, smaller median time delays for 
small opening angle pairs, larger median time delays for 
asymmetric pairs, can be seen in this Figure as well. In 
addition, the scatter presented here is useful to check 
which lens systems are less dependent on lens potentials: 
Our result indicates that asymmetric image pairs that 
are collinear with their lens galaxies (6*^ ~ 180°) have 
the least scatter and therefore are more suitable for the 
determination of the Hubble constant. It is interesting to 
note that in the statistical sense double lenses are more 
valuable than quadruple lenses because it shows a smaller 
sensitivity to various lens potentials. In addition, double 
lenses have advantages of easier measurements of time 
delays in observations and smaller fractional errors (be- 
cause double lenses have longer time delays than quads). 
This is in contrast to individual mass modeling in which 
quadruple lenses are more useful because of much more 
observational constraints on mass models. 

5. COMPARISON WITH OBSERVED TIME DELAYS 

We now examine if the conditional distribution com- 
puted in ^ is consistent with the observed distribution 
of time delays. Since we need to assume the Hubble con- 
stant in order to convert observed time delays to reduce 
time delays S, in this section we adopt h — 0.73 that was 
obtained from the combined analysis of the cosmic mi- 
crowave background (CMB) a nisotropy and clustering of 
galaxies (jTegmark et al.|[2"006[ ). We derive reduced time 
delays for the 17 published time delay quasars (41 image 
pairs), which is summarized in Table [T] 

First, we compare these reduced time delays with the 
conditional probability plotted in Figure [H We show 
p{E\Rij) and p{E\9ij) with observed values overplotted 
in Figure [S) It appears that the data are roughly consis- 
tent with our probability distribution from simulations. 
The large scatter of very symmetric lenses is also seen 
for observed time delays. It is interesting that observed 
time delays appear to exhibit the decline of time delays at 



small opening angles, just as our theoretical model pre- 
dicts. We may need more quasar time delays to confirm 
this trend observationally. 

More directly, observed time delays should be com- 
pared with the probability distribution for given Rij and 
9ij, i.e., p{^\Rij,6ij). In Figure[7]we compare observed 
S with the probability distribution of S expected from 
corresponding image configuration of each image pair. 
We find that the distributions agree with observed val- 
ues on average, similarly as Figure [SI The reduced time 
delays of B0218-h357 and SBS0909+532 have larger er- 
rors because the positions of the lens galaxies are poorly 
determined. The large error of HE0435— 1223 AC comes 
from small — rl, whereas that of B 1422-1-231 is sim- 
ply because of the large time delay measurement errors. 
We comment that RXJ0911-^0551 and Q0957-H561 has 
significantly smaller S compared with our model predic- 
tion. This is clearly because of the cluster convergence: 
The two lenses lie in near the centers of clusters and 
therefore the observed time delays are pushed down by 
the convergence coming from dark matter in the clusters 
(see Sj6]). The reduced time delays of B1608-I-656 and 
SBS1520-I-530 are largely offset from the predicted val- 
ues, and this is probably because of satellite galaxies in 
the lens systems that significantly affect the time delays. 
The large tim e delays of RXJ1131- 1231 AB and AC were 
also noted by [Morgan et al.l (|2007D : Our result indicates 
that the broadened theoretical distributions due to small 
perturbations are enough to explain the observed high 
values of the time delays. The large offsets of B1422-I-231 
and SDSS J1650-H4251 may come from the large uncer- 
tainties of measured time delays (see ^|6]). 

6. IMPLICATIONS FOR Ho 

In this section, we turn the problem around and con- 
strain the Hubble constant using the conditional prob- 
ability distribution constructed from the Monte-Carlo 
simulations. Although the current sample of observed 
time delays summarized in Table [T] is somewhat hetero- 
geneous and may not be appropriate for the statistical 
study, we do this to demonstrate how we can constrain 
the Hubble constant from the probability distribution. 
We basically take all lens systems listed in Table [l] but 
adopt the following setup to reduce systematic errors. 

• We do not use the AB time delay of SDSS 
J1004-I-4112. It is a lens system caused by a mas- 
sive cluster of galaxies, thus our input distribution, 
which is designed for galaxies, does not represent a 
fair distribution of lens potentials. Moreover, the 
center of the lens potential appears to be offset 
from the position of the brightest cluster galaxy 
([Oguri et al.ll2004j ). which makes it inaccurate to 
estimate important parameters such S, Rij, and 
Bij. 

• Both RXJ091H-G551 and QG957-f 561 are known to 
reside in the cluster environment, thus they are sig- 
nificantly affected by the cluster convergence Kck,- 
Since the mass-sheet degeneracy says S c>c At (x 
1 — Kciu, we divide reduce time delays S for these 
two systems by 1 — KcIu in order to deconvolve the 
effect of the cluster convergence. As the values of 
Kciu, we adopt Kdu = 0.3±0.04 for RXJ091 1-1-0551 
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Summary of Observed Quasar Time Delays 
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[2002; (12) CASTLES (littp://cfa-www.harvard.cdu/castlcs/ ); (13) Oscoz ct al. 1997; (14) Lubin ct al. 2000; (15) Ullan et al. 20Q|; 
(16) Schcchtcr ct al. 1998; (17) Jakobsson et al. 2005; (18) Eigcnbrod ct al. 2007; (19) Walsh ct al. 1979; (20) Young et al. 1981; (21) 
[Kundic et al. 1997; (22) Barkana et al. 1999; (23) Inada et al. 2003; (24) Oguri et al. 2004; (25) Inada ct al. 2005; (26) Fohlmcistcr' et aU 
[2007; (27) Wisotzki ct al. 1993; (28) Lidman ct al. 2000; (29) Ofek & Maoz 2003; (30) Poindcxter et al. 2006; (31) Wcymann et al. lOStt 
(32) Schcchtcr et al. 1997; (33) Barkana 1997; (34) Tonrv 1998; (35) Impcy ct al. 1998; (36) Chartas et al. 2004; (37) Sluse ct al. 200l; 
(38) Morgan ct al. 2007; (39) Patnaik et al. 1992; (40) Patnaik & Narasimha 2001; (41) Cliavusliyan et al. 1997; (42) Burud et al. 20022; 
j43) Faurc et al. 2002; (44) Jackson et al. 1995; (45) Fassnaelit & Cohen 1998; (46) Koopmans ct al. 1998; (47) Burud ct al. 2000; (48) 
[Myers ct al. 1995; (49) Fassnacht ct al. 1996; (50) Koopmans & Fassnacht 1999; (51) Fassnacht et al. 2002; (52) Koopmans et al. 2003; (53) 
^rgan et al..,2003i : (54) iVuissoz et al.,,2007.: (55) igubrahmanvan et al., .199Q: (56) .Wiklind fc Combes.. 1996 : (57) .Lovell et al..,1998i ; (58) 
iLidman et anil99g ; (59) IWisotzki et'aLlll996l ; (60) [Eirud et al. 20024 

Note. — All measured time delays are listed, except time delays between images Al — A3 of RXJ091 1+0551 and images A — D of 
Q2237+030 (Vakulik ct al. 2006) for which error-bars are large and there fore the de t ection s are marginal. For Q2237+030, a possible X-ray 
detection of the time delay between image A and B was also reported by lDai et al.l II2003I ). Errors indicate Itr. 

The values of reduced time delays H computed from observed time delays and image configurations. For the Hubble constant h = 0.73 is 
assumed. Errors of S come from those of Ai and — ffl- ^ We assume the position of the brightest cluster galaxy Gl for the center of 
the lens potential, though mass modeling implies the significant offset of the potential center from Gl l lOguri et al.ll2003) .'^ The values and 
errors of time delays were not directly given in the literature, thus we inferred them from those of other image pairs. 
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Fig. 7. — The conditional probability distributions p(H|i?ij, 6ij) are compared with reduced time delays H computed from observed time 
delays of 17 published time delay quasars. Again, the Hubble constant o{ h = 0.73 is assumed. Solid line curves are in fact differential 
probability distribution . In each panel, we draw the differential probability distribution dp/dlogE, computed from observed values of Rij 
and 8ij by a solid line. Reduced time delays from observed time delays are shown by vertical dashed lines plus error- bars. See Table [Tl for 
the summary of observations. 



Gravitational Lens Time Delays 



13 




Fig. 8. — Statistical constraint on the Hubble constant from 
16 time delay quasars (40 image pairs). Thick solid line indicates 
goodness-of-fit parameter from all 16 lens systems plotted as a 
function of the Hubble constant h. The resulting Hubble constant 
is /i = 0.70^ ^ "'^ 



^ 0.70to 05 at 95% confi- 



-'io 02 68% confidence and h - 
dence. The Hubble constant estimated using jackknife resampling 
has a larger error, h = 0.70 ± 0.06 at 68% confidence (see text for 
details). Thin solid lines show goodness-of-fit parameter for each 
lens system. 



HHiorth et al l [200^. and kh., = 0.26 ± 0.08 for 
QG957-H561 (iFischer et al.lfl997l) . 

In summary, we use 16 lensed quasar systems (40 image 
pairs) to constrain the Hubble constant. For each image 
pair, we compute the likelihood as follows: 



Cpih) 



/dp 
— (S|i?ij^obs, Sij,obs)G(E\Eohs{h))dE, 



(25) 



where i?y^obs, %,obsi Sobs are those for this specific image 
pair listed in Table [H and G'(S|E!obs('i)) indicates the 
Gaussian distribution with median S = Sobs- Note that 
calculating Sobs from observed time delays require the 
Hubble constant h, hence Cp is a function of h. Then 
we compute the effective chi-square by summing up the 
logarithm of the likelihoods: 



E ;^E[-21n/:.(/^)] 



(26) 



quasar 



The first summation runs over lens systems, whereas the 
second summation runs over image pairs for each lens 
system; the number of pairs for each lens is denoted by 
Up. Note that Up = 1 all double lens systems, and a 
quadruple lens system should have np < 4C2 = 6 depend- 
ing on how many time delays have been observed for the 
lens system. The factor l/np was introduced such that all 
lens systems have equal weight on the effective chi-square 
irrespective of the number of image pairs. We derive the 
best-fit value and its error of h by the standard way using 
a goodness-of-fit parameter Axcff = Xcff — Xoff(inin). 

We show our result in Figure [H The Hubble constant 
measured from the combination of all 16 lens systems is 



0.70] 



0.70 



+0.09 



at 



^^0 ^02 at 68% confidence and h — u.iu_q o5 
95% confidence. The obtained value is in good agree- 
ment with other estimates, such as th e local distance 
meas urement using Ccpheid calibr ation (iFreedman et al.l 
I2OOID and the CMB anisotropy (jTegmark et all 120061 : 



TABLE 2 
Hubble Constant from Each Lens 
System 



Lens Name 



h (la range) 



B0218-I-357 

HE0435-1223 

RXJ0911+0551 

SBS0909-f532 

FBQ095H-2635 

Q0957+561 

HE1104-1805 

PG1115-I-080 

RXJ1131-1231 

B1422-I-231 

SBS1520-I-530 

B1600-(-434 

B1608-(-656 

SDSS J1650+4251 

PKS1830-211 

HE2149-2745 



0.21 (-) 
1.02 (0.70-1.39) 
0.96 (0.75-1.21) 

0.84 (0.47-) 
0.67 (0.56-0.81) 
0.99 (0.82-1.17) 
1.04 (0.92-1.22) 
0.66 (0.49-0.84) 
0.79 (0.59-1.03) 

0.16 (-0.36) 
0.53 (0.46-0.61) 
0.65 (0.54-0.77) 
0.89 (0.77-1.20) 
0.53 (0.44-0.63) 

0.88 (0.58-) 
0.69 (0.57-0.82) 



0.70 (0.68-aW 



Note. — The Hubble constant and 
its error are estimated from the effec- 
tive chi-square. 



ISpergel et al.ll2007f l. The constraint from each lens sys- 
tem, which is plotted in Figure[51 is summarized in Table 

m 

We also derive the Hubble constant using the jack- 
knife resampling by regarding each 16 lens system as 
a subsample. The result h = 0.70 ± 0.06 at 68% con- 
fidence has the same mean but larger error than that 
estimated from the effective chi-square. There are sev- 
eral possible source of this difference. One is the un- 
derestimate of the width of the input distributions. In 
particular, many of the time delay quasar systems has 
been claime d to be affected by len s galaxy environ- 
ments (e.g.. [Morgan et al.l 120051: iFass nacht ct al. 200^ 



Momcheva et al. I I2006I: IWiUiams et al.ll2006l : I Auger et aU 



2007D . and thus our input strength of external shear 
might be somewhat smaller than the true one (see also 
discussion in fjT]). Another possible source is the non- 
Gaussianity of measured time delays: In equation (j25p 
we assumed the Gaussian distribution for the measure- 
ment uncertainties of time delays, but sometimes they 
are quite different from the Gaussian distribution."* We 
note that in our method we can in principle include non- 
Gaussianity by just replacing G(S) in equation (j25p with 
any appropriate probability distributions, as long as we 
know such distributions. 

7. DISCUSSIONS AND CONCLUSION 

In this paper, we have studied time delays between 
multiply imaged quasars. Adopting the reduced time de- 
lay, which is a measure of how the lens potential is com- 
plicated compared with the simple isothermal form, we 
have explored the dependence of time delays on various 
complex structure of lens potentials such as external per- 
turbations, non-isothermality, and substructures. The 
distribution of time delays has been studied as a func- 

Among time delays listed Table [T] those of SDSS J1650-(-4251 
and B1422+231 could be significantly different from the true values 
(C. S. Kochanek, private communication). We perform the same 
analysis excluding these two systems and find the Hubble constant 
to be h = 0.70io o4 at 68% confidence from the effective chi-square. 
Therefore our result is not biased significantly by these systems. 
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tion of image configuration which we characterize using 
two dimensionless quantities, the asymmetry and open- 
ing angle of an image pair. We have pointed out that the 
sensitivity on lens potentials is quite dependent on the 
image configuration. For instance, more symmetric im- 
age pairs are more affected by a small change of the lens 
potential. Image pairs with smaller opening angles are 
also more sensitive to lens potentials. In particular time 
delays of close image pairs are very sensitive to higher- 
order external perturbations and substructures that are 
very hard to be constrained from mass modeling even for 
best studied lens systems. In addition, those pairs usu- 
ally have the largest relative uncertainties of time delay 
measurements. Therefore we conclude that it is quite 
difficult to extract any information on the Hubble con- 
stant or the mass distribution from close image pairs. It 
is interesting to note that perturbations on lens potential 
not only introduce scatter around the mean but also can 
systematically shift the distribution of time delays. One 
such example is smaller time delays for smaller opening 
angle image pairs, which is caused by external shear. 

We have performed Monte-Carlo simulations to de- 
rive a probability distribution of reduced time delays 
for each image configuration. Input distributions are de- 
termined from observational and theoretical constraints. 
The distribution is weighted by the lensing cross sec- 
tion and magnification bias, allowing a realistic estimate 
of time delay distributions. The probability distribu- 
tion was then compared with observed time delays. We 
have shown that the distribution of time delays com- 
puted from our simulations is in good agreement with 
observed time delays. In particular, distributions of ob- 
served time delays also exhibit strong dependence of im- 
age configuration in a consistent manner with our theo- 
retical expectations. The probability distribution can be 
used to constrain the Hubble constant. We have found 
that 16 published time delay quasars constrain it to be 
h = 0.70j^Q;Q2 at 68% confidence using the effective chi- 
square or h = 0.70 ± 0.06 estimated using jackknife re- 
sampling, consistent with other estimates. 

An important caveat is that our lensed quasar sam- 
ple is quite heterogeneous. In particular, it should be 
noted that current time delay quasars (see Table [T]) have 
significantly larger image separations on average com- 
pared with the other quasar lenses: The median im- 
age separation of time delay quasars listed in Table [1] 
is 1'.'7 (image separations before and after the median 
are 1'.'5 and 2'.'1), whereas that of all lensed quasars 
is ~ 1'.'4. Quasar lenses with larger image separa- 
tions are more likely to lie in dense environments be- 
cause both the image separation and biased cross section 
are boosted by su rrounding dark matter (Ogur i et al.l 
I2005t iQguril |2006| ). thus the Hubble constant inferred 
from those lenses are more affected by the environmen- 
tal convergence. Indeed, the association of group/cluster 
has been rep orted for mor e thaii half of the time delay 
quasars (e.g., [M organ ct al. 2005; Fassnacht ct al. 2006; 



the effective chi-square is consistent with our full result, 
thus we conclude that the effect of lens galaxy environ- 
ments is not so drastic here. However, to minimize the 
systematic effect, in the future we should apply our sta- 
tistical method to well-defin ed samples of lensed quasars 
such as the CLASS (Mvers et all [200l iBrowne et all 



Moincheva et al. 2006 : Williams c^aL 2006; Aug er ct al.. 
20071) . In addition, our input distributions of external 



perturbations may be underestimated for these large sep- 
aration lenses. To estimate this systematic effect, we 
exclude four lensed quasars with image separation larger 
than 3" and repeat the analysis done in fj6] The resulting 
Hubble constant h = 0.67tn at 68% confidence from 



I2003D . SQLS (jOguri et al.ll2006D . and those obtained in 
future lens surveys. 

Another source of the systematic effect is the un- 
certainty of our input distribution of lens potentials. 
Among others, the most important systematic error 
comes from the uncertainty of the mean value of the 
slope of the radial profile, a. While a = 1 (isother- 
mal) for the mean appears to be a reasonable choice, 
direc t studies of len s galaxies [e.g.. iTreu fc KoopmansI 
12004 iRusin fc Ko cha nek 20051: iHamana et al.ll2007| ) in- 
dicate that the error on the mean a could be as large as 
0.1. The derived Hubble constant depends on the slope 
as h oc 2 — a, therefore the change of the mean a sys- 
tematically shifts the best-fit value. The scaling relation 
suggests that the 0.1 error of the mean results in 10% 
error on h, indicating that the systematic error may be 
even larger than the statistical error. Another important 
systematic effect is caused by external convergence from 
lens galaxy environments. Since its effect o n the Hub- 
ble co nstant is straightforward (e.g., Keeton fc Zabludo3 
l2004f ). one can estimate the effect of external convergence 
rather easily even without includin g the distribution i n 
the simulation. From the result of I Oguri et al.l (|2005( l. 
we expect that the posterior distribution of external con- 
vergence with lensing bias taken into account is roughly 
Kcxt ~ 0.03±0.03 unless the image separation is too large. 
This, combined with the fact that the Hubble constant 
scales as h (X 1 — K^xt , suggest that the effect of external 
convergence is not so dominant here (note that exter- 
nal convergence was already taken into account for two 
extreme lenses, Q0957-I-561 and RXJ091H~0551). By in- 
cluding rough estimates of these two systematic errors, 
we obtain h = 0.68 ± 0.06(stat.) ± 0.08(syst.), indicating 
the importance of reducing the systematic error. 

Since the Hubble constant is now determined fairly 
well by other methods, time delays are sometimes used 
to study mass distributions of lens objects. Our statis- 
tical technique offers a new method to study the lens 
mass distribution. By comparing the probability distri- 
butions for different input distributions of lens potentials 
(e.g., different median slopes of the primary lens galaxy), 
one can infer which input model is most plausible. Un - 
like previous statistical studies (e.g.. lOguri et al.ll2002f) . 
this new method allows us to include various complex- 
ity of lens potentials relatively easily, particularly the 
non-spherically symmetric nature of lens potentials. We 
note that our input distribution of lens potentials was 
designed for galaxies, but it is straightforward to modify 
it to that of lensing by other populations, e.g., lensing by 
a cluster of galaxies. 

In summary, our new statistical approach is invaluable 
for the study of both cosmological parameters and struc- 
ture of lens potentials. We believe its importance grows 
more and more in the era of large-scale synoptic surveys 
such as LSST and SNAP: Quasar lens candidates are 
easily recognized in these synoptic surve ys by rnaking 
use of strong time var iability of quasars (jPindorl l2005l : 
iKochanek et al.ll2006a^ . Strong lensing of distant super- 
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novae offers additional interesting opportunity to apply 
our statistical technique. The statistical analysis is es- 
sential to make efficient use of the large homogeneous 
samples of strong lenses provided by these surveys. 
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cussions and comments. I am grateful to an anonymous 
referee for his/her many useful suggestions. This work 
was supported in part by the Department of Energy con- 
tract DE-AC02-76SF00515. 
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